####################################
# Code to Replicate Charts Figures and Tables
# Trust in Public Policy Algorithms
# Journal of Politics
# Criminal Justice Experiments
# Study 2
# Ryan Kennedy, Philip D. Waggoner, and Matthew Ward
# 2021
####################################

#### Load Needed Libraries
library(arm)
library(gridExtra)
library(tidyverse)
library(here)
library(stargazer)
library(gridExtra)


#### Load data
cj1 <- read_csv(here("SentencingExperiment1.csv"))
cj2 <- read_csv(here("SentencingExperiment2.csv"))

#### Create variables and controls
# Add in the information on advice and control variables
cj1 <- cj1 %>%
  mutate(adviceSC2 = 80, adviceSC4 = 68, adviceSC7 = 52, adviceSC11 = 32,
         adviceSC12 = 56, adviceSC13 = 36, adviceSC18 = 37, adviceSC19 = 87,
         advicebSC2 = 100 - 80, advicebSC4 = 100 - 68, advicebSC7 = 100 - 52, advicebSC11 = 100 - 32,
         advicebSC12 = 100 - 56, advicebSC13 = 100 - 36, advicebSC18 = 100 - 37, advicebSC19 = 100 - 87,
         encouragement = ifelse(!is.na(treatment_explain), 1, 0),
         female = ifelse(gender == "Female", 1, 0),
         edvalue = case_when(education == "Elementary or some high school" ~ 1,
                             education == "High school graduate/GED" ~ 2,
                             education == "Trade or vocational certification" ~ 2,
                             education == "Some college, or an associate degree" ~ 3,
                             education == "Bachelor's degree (for example, BA, AB, BS)" ~ 3,
                             education == "Some graduate school" ~ 4,
                             education == "Master's degree (for example, MA, MSW, MBA)" ~ 4,
                             education == "Professional degree (for example, MD, JD, DDS)" ~ 4,
                             education == "Doctoral degree (for example, PhD, EdD)" ~ 4),
         toa1 = case_when(trustinautomation_1 == "Strongly agree" ~ 7,
                          trustinautomation_1 == "Agree" ~ 6,
                          trustinautomation_1 == "Somewhat agree" ~ 5,
                          trustinautomation_1 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_1 == "Somewhat disagree" ~ 3,
                          trustinautomation_1 == "Disagree" ~ 2,
                          trustinautomation_1 == "Strongly disagree" ~ 1),
         toa2 = case_when(trustinautomation_2 == "Strongly agree" ~ 7,
                          trustinautomation_2 == "Agree" ~ 6,
                          trustinautomation_2 == "Somewhat agree" ~ 5,
                          trustinautomation_2 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_2 == "Somewhat disagree" ~ 3,
                          trustinautomation_2 == "Disagree" ~ 2,
                          trustinautomation_2 == "Strongly disagree" ~ 1),
         toa3 = case_when(trustinautomation_3 == "Strongly agree" ~ 7,
                          trustinautomation_3 == "Agree" ~ 6,
                          trustinautomation_3 == "Somewhat agree" ~ 5,
                          trustinautomation_3 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_3 == "Somewhat disagree" ~ 3,
                          trustinautomation_3 == "Disagree" ~ 2,
                          trustinautomation_3 == "Strongly disagree" ~ 1),
         toa4 = case_when(trustinautomation_4 == "Strongly agree" ~ 1,
                          trustinautomation_4 == "Agree" ~ 2,
                          trustinautomation_4 == "Somewhat agree" ~ 3,
                          trustinautomation_4 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_4 == "Somewhat disagree" ~ 5,
                          trustinautomation_4 == "Disagree" ~ 6,
                          trustinautomation_4 == "Strongly disagree" ~ 7),
         toa5 = case_when(trustinautomation_5 == "Strongly agree" ~ 1,
                          trustinautomation_5 == "Agree" ~ 2,
                          trustinautomation_5 == "Somewhat agree" ~ 3,
                          trustinautomation_5 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_5 == "Somewhat disagree" ~ 5,
                          trustinautomation_5 == "Disagree" ~ 6,
                          trustinautomation_5 == "Strongly disagree" ~ 7),
         toa6 = case_when(trustinautomation_6 == "Strongly agree" ~ 7,
                          trustinautomation_6 == "Agree" ~ 6,
                          trustinautomation_6 == "Somewhat agree" ~ 5,
                          trustinautomation_6 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_6 == "Somewhat disagree" ~ 3,
                          trustinautomation_6 == "Disagree" ~ 2,
                          trustinautomation_6 == "Strongly disagree" ~ 1),
         toa7 = case_when(trustinautomation_7 == "Strongly agree" ~ 7,
                          trustinautomation_7 == "Agree" ~ 6,
                          trustinautomation_7 == "Somewhat agree" ~ 5,
                          trustinautomation_7 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_7 == "Somewhat disagree" ~ 3,
                          trustinautomation_7 == "Disagree" ~ 2,
                          trustinautomation_7 == "Strongly disagree" ~ 1),
         toa8 = case_when(trustinautomation_8 == "Strongly agree" ~ 1,
                          trustinautomation_8 == "Agree" ~ 2,
                          trustinautomation_8 == "Somewhat agree" ~ 3,
                          trustinautomation_8 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_8 == "Somewhat disagree" ~ 5,
                          trustinautomation_8 == "Disagree" ~ 6,
                          trustinautomation_8 == "Strongly disagree" ~ 7),
         trustAutomation = (toa1 + toa2 + toa3 + toa4 + toa5 + toa6 + toa7 + toa8)/8,
         partisanship = case_when(partylean == "I lean neither way" ~ 4,
                                  partylean == "Lean Democrat" ~ 3,
                                  partylean == "Lean Republican" ~ 5,
                                  party == "Republican" & partystrong == "No, NOT strongly" ~ 6,
                                  party == "Democrat" & partystrong == "No, NOT strongly" ~ 2,
                                  party == "Republican" & partystrong == "Yes, strongly" ~ 7,
                                  party == "Democrat" & partystrong == "Yes, strongly" ~ 1),
         needcog1 = case_when(needforcognition_1 == "Extremely characteristic" ~ 5,
                              needforcognition_1 == "Somewhat characteristic" ~ 4,
                              needforcognition_1 == "Uncertain" ~ 3,
                              needforcognition_1 == "Somewhat uncharacteristic" ~ 2,
                              needforcognition_1 == "Extremely uncharacteristic" ~ 1),
         needcog2 = case_when(needforcognition_2 == "Extremely characteristic" ~ 5,
                              needforcognition_2 == "Somewhat characteristic" ~ 4,
                              needforcognition_2 == "Uncertain" ~ 3,
                              needforcognition_2 == "Somewhat uncharacteristic" ~ 2,
                              needforcognition_2 == "Extremely uncharacteristic" ~ 1),
         needcog3 = case_when(needforcognition_3 == "Extremely characteristic" ~ 5,
                              needforcognition_3 == "Somewhat characteristic" ~ 4,
                              needforcognition_3 == "Uncertain" ~ 3,
                              needforcognition_3 == "Somewhat uncharacteristic" ~ 2,
                              needforcognition_3 == "Extremely uncharacteristic" ~ 1),
         needcog4 = case_when(needforcognition_4 == "Extremely characteristic" ~ 1,
                              needforcognition_4 == "Somewhat characteristic" ~ 2,
                              needforcognition_4 == "Uncertain" ~ 3,
                              needforcognition_4 == "Somewhat uncharacteristic" ~ 4,
                              needforcognition_4 == "Extremely uncharacteristic" ~ 5),
         needjudge = (needcog1 + needcog3)/2,
         needcognition = (needcog2 + needcog4)/2,
         source = "MTurk")

cj2 <- cj2 %>%
  mutate(adviceSC2 = 80, adviceSC4 = 68, adviceSC7 = 52, adviceSC11 = 32,
         adviceSC12 = 56, adviceSC13 = 36, adviceSC18 = 37, adviceSC19 = 87,
         advicebSC2 = 100 - 80, advicebSC4 = 100 - 68, advicebSC7 = 100 - 52, advicebSC11 = 100 - 32,
         advicebSC12 = 100 - 56, advicebSC13 = 100 - 36, advicebSC18 = 100 - 37, advicebSC19 = 100 - 87,
         encouragement = ifelse(!is.na(treatment_explain_1), 1, 0),
         female = ifelse(gender == 2, 1, 0),
         edvalue = case_when(education == 1 ~ 1,
                             education == 2 ~ 2,
                             education == 3 ~ 2,
                             education == 4 ~ 3,
                             education == 5 ~ 3,
                             education == 6 ~ 3,
                             education == 7 ~ 4,
                             education == 8 ~ 4),
         partisanship = case_when(political_party == 7 | political_party == 4 ~ 4,
                                  political_party == 6 | political_party == 3 ~ 3,
                                  political_party == 8 | political_party == 5 ~ 5,
                                  political_party == 9 ~ 6,
                                  political_party == 2 ~ 2,
                                  political_party == 10 ~ 7,
                                  political_party == 1 ~ 1),
         toa1 = case_when(trustinautomation_1 == "Strongly agree" ~ 7,
                          trustinautomation_1 == "Agree" ~ 6,
                          trustinautomation_1 == "Somewhat agree" ~ 5,
                          trustinautomation_1 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_1 == "Somewhat disagree" ~ 3,
                          trustinautomation_1 == "Disagree" ~ 2,
                          trustinautomation_1 == "Strongly disagree" ~ 1),
         toa2 = case_when(trustinautomation_2 == "Strongly agree" ~ 7,
                          trustinautomation_2 == "Agree" ~ 6,
                          trustinautomation_2 == "Somewhat agree" ~ 5,
                          trustinautomation_2 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_2 == "Somewhat disagree" ~ 3,
                          trustinautomation_2 == "Disagree" ~ 2,
                          trustinautomation_2 == "Strongly disagree" ~ 1),
         toa3 = case_when(trustinautomation_3 == "Strongly agree" ~ 7,
                          trustinautomation_3 == "Agree" ~ 6,
                          trustinautomation_3 == "Somewhat agree" ~ 5,
                          trustinautomation_3 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_3 == "Somewhat disagree" ~ 3,
                          trustinautomation_3 == "Disagree" ~ 2,
                          trustinautomation_3 == "Strongly disagree" ~ 1),
         toa4 = case_when(trustinautomation_4 == "Strongly agree" ~ 1,
                          trustinautomation_4 == "Agree" ~ 2,
                          trustinautomation_4 == "Somewhat agree" ~ 3,
                          trustinautomation_4 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_4 == "Somewhat disagree" ~ 5,
                          trustinautomation_4 == "Disagree" ~ 6,
                          trustinautomation_4 == "Strongly disagree" ~ 7),
         toa5 = case_when(trustinautomation_5 == "Strongly agree" ~ 1,
                          trustinautomation_5 == "Agree" ~ 2,
                          trustinautomation_5 == "Somewhat agree" ~ 3,
                          trustinautomation_5 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_5 == "Somewhat disagree" ~ 5,
                          trustinautomation_5 == "Disagree" ~ 6,
                          trustinautomation_5 == "Strongly disagree" ~ 7),
         toa6 = case_when(trustinautomation_6 == "Strongly agree" ~ 7,
                          trustinautomation_6 == "Agree" ~ 6,
                          trustinautomation_6 == "Somewhat agree" ~ 5,
                          trustinautomation_6 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_6 == "Somewhat disagree" ~ 3,
                          trustinautomation_6 == "Disagree" ~ 2,
                          trustinautomation_6 == "Strongly disagree" ~ 1),
         toa7 = case_when(trustinautomation_7 == "Strongly agree" ~ 7,
                          trustinautomation_7 == "Agree" ~ 6,
                          trustinautomation_7 == "Somewhat agree" ~ 5,
                          trustinautomation_7 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_7 == "Somewhat disagree" ~ 3,
                          trustinautomation_7 == "Disagree" ~ 2,
                          trustinautomation_7 == "Strongly disagree" ~ 1),
         toa8 = case_when(trustinautomation_8 == "Strongly agree" ~ 1,
                          trustinautomation_8 == "Agree" ~ 2,
                          trustinautomation_8 == "Somewhat agree" ~ 3,
                          trustinautomation_8 == "Neither agree nor disagree" ~ 4,
                          trustinautomation_8 == "Somewhat disagree" ~ 5,
                          trustinautomation_8 == "Disagree" ~ 6,
                          trustinautomation_8 == "Strongly disagree" ~ 7),
         trustAutomation = (toa1 + toa2 + toa3 + toa4 + toa5 + toa6 + toa7 + toa8)/8, 
         needcog1 = case_when(needforcognition_1 == "Extremely characteristic" ~ 5,
                              needforcognition_1 == "Somewhat characteristic" ~ 4,
                              needforcognition_1 == "Uncertain" ~ 3,
                              needforcognition_1 == "Somewhat uncharacteristic" ~ 2,
                              needforcognition_1 == "Extremely uncharacteristic" ~ 1),
         needcog2 = case_when(needforcognition_2 == "Extremely characteristic" ~ 5,
                              needforcognition_2 == "Somewhat characteristic" ~ 4,
                              needforcognition_2 == "Uncertain" ~ 3,
                              needforcognition_2 == "Somewhat uncharacteristic" ~ 2,
                              needforcognition_2 == "Extremely uncharacteristic" ~ 1),
         needcog3 = case_when(needforcognition_3 == "Extremely characteristic" ~ 5,
                              needforcognition_3 == "Somewhat characteristic" ~ 4,
                              needforcognition_3 == "Uncertain" ~ 3,
                              needforcognition_3 == "Somewhat uncharacteristic" ~ 2,
                              needforcognition_3 == "Extremely uncharacteristic" ~ 1),
         needcog4 = case_when(needforcognition_4 == "Extremely characteristic" ~ 1,
                              needforcognition_4 == "Somewhat characteristic" ~ 2,
                              needforcognition_4 == "Uncertain" ~ 3,
                              needforcognition_4 == "Somewhat uncharacteristic" ~ 4,
                              needforcognition_4 == "Extremely uncharacteristic" ~ 5),
         needjudge = (needcog1 + needcog3)/2,
         needcognition = (needcog2 + needcog4)/2,
         source = "Lucid")

# Record the source of the advice respondent received
cj1 <- cj1 %>%
  mutate(adTypeSC2 = ifelse(!is.na(sc2a_1) | !is.na(s2aoa_1), "algorithm", 
                            ifelse(!is.na(sc2j_1) | !is.na(s2joa_1), "judge",
                                   ifelse(!is.na(sc2mt_1) | !is.na(s2moa_1), "mt", NA))),
         adTypeSC4 = ifelse(!is.na(sc4a_1) | !is.na(s4aoa_1), "algorithm", 
                            ifelse(!is.na(sc4j_1) | !is.na(s4joa_1), "judge",
                                   ifelse(!is.na(sc4mt_1) | !is.na(s4moa_1), "mt", NA))),
         adTypeSC7 = ifelse(!is.na(sc7a_1) | !is.na(s7aoa_1), "algorithm", 
                            ifelse(!is.na(sc7j_1) | !is.na(s7joa_1), "judge",
                                   ifelse(!is.na(sc7mt_1) | !is.na(s7moa_1), "mt", NA))),
         adTypeSC11 = ifelse(!is.na(sc11a_1) | !is.na(s11aoa_1), "algorithm", 
                             ifelse(!is.na(sc11j_1) | !is.na(s11joa_1), "judge",
                                    ifelse(!is.na(sc11mt_1) | !is.na(s11moa_1), "mt", NA))),
         adTypeSC12 = ifelse(!is.na(sc12a_1) | !is.na(s12aoa_1), "algorithm", 
                             ifelse(!is.na(sc12j_1) | !is.na(s12joa_1), "judge",
                                    ifelse(!is.na(sc12mt_1) | !is.na(s12moa_1), "mt", NA))),
         adTypeSC13 = ifelse(!is.na(sc13a_1) | !is.na(s13aoa_1), "algorithm", 
                             ifelse(!is.na(sc13j_1) | !is.na(s13joa_1), "judge",
                                    ifelse(!is.na(sc13mt_1) | !is.na(s13moa_1), "mt", NA))),
         adTypeSC18 = ifelse(!is.na(sc18a_1) | !is.na(s18aoa_1), "algorithm", 
                             ifelse(!is.na(sc18j_1) | !is.na(s18joa_1), "judge",
                                    ifelse(!is.na(sc18mt_1) | !is.na(s18moa_1), "mt", NA))),
         adTypeSC19 = ifelse(!is.na(sc19a_1) | !is.na(s19aoa_1), "algorithm", 
                             ifelse(!is.na(sc19j_1) | !is.na(s19joa_1), "judge",
                                    ifelse(!is.na(sc19mt_1) | !is.na(s19moa_1), "mt", NA))))

cj2 <- cj2 %>%
  mutate(adTypeSC2 = ifelse(!is.na(s2aoa_1), "algorithm", 
                            ifelse(!is.na(s2joa_1), "judge",
                                   ifelse(!is.na(s2moa_1), "mt", NA))),
         adTypeSC4 = ifelse(!is.na(s4aoa_1), "algorithm", 
                            ifelse(!is.na(s4joa_1), "judge",
                                   ifelse(!is.na(s4moa_1), "mt", NA))),
         adTypeSC7 = ifelse(!is.na(s7aoa_1), "algorithm", 
                            ifelse(!is.na(s7joa_1), "judge",
                                   ifelse(!is.na(s7moa_1), "mt", NA))),
         adTypeSC11 = ifelse(!is.na(s11aoa_1), "algorithm", 
                             ifelse(!is.na(s11joa_1), "judge",
                                    ifelse(!is.na(s11moa_1), "mt", NA))),
         adTypeSC12 = ifelse(!is.na(s12aoa_1), "algorithm", 
                             ifelse(!is.na(s12joa_1), "judge",
                                    ifelse(!is.na(s12moa_1), "mt", NA))),
         adTypeSC13 = ifelse(!is.na(s13aoa_1), "algorithm", 
                             ifelse(!is.na(s13joa_1), "judge",
                                    ifelse(!is.na(s13moa_1), "mt", NA))),
         adTypeSC18 = ifelse(!is.na(s18aoa_1), "algorithm", 
                             ifelse(!is.na(s18joa_1), "judge",
                                    ifelse(!is.na(s18moa_1), "mt", NA))),
         adTypeSC19 = ifelse(!is.na(s19aoa_1), "algorithm", 
                             ifelse(!is.na(s19joa_1), "judge",
                                    ifelse(!is.na(s19moa_1), "mt", NA))))

# Record the predictions made under each scenario by the participant
# Note: in this case, the rowMeans() only records the value for the one they filled in, since all others will be NA
cj1 <- cj1 %>%
  mutate(predSC2 = rowMeans(cbind(sc2a_1, s2aoa_1, sc2j_1, s2joa_1, sc2mt_1, s2moa_1), na.rm = TRUE),
         predSC4 = rowMeans(cbind(sc4a_1, s4aoa_1, sc4j_1, s4joa_1, sc4mt_1, s4moa_1), na.rm = TRUE),
         predSC7 = rowMeans(cbind(sc7a_1, s7aoa_1, sc7j_1, s7joa_1, sc7mt_1, s7moa_1), na.rm = TRUE),
         predSC11 = rowMeans(cbind(sc11a_1, s11aoa_1, sc11j_1, s11joa_1, sc11mt_1, s11moa_1), na.rm = TRUE),
         predSC12 = rowMeans(cbind(sc12a_1, s12aoa_1, sc12j_1, s12joa_1, sc12mt_1, s12moa_1), na.rm = TRUE),
         predSC13 = rowMeans(cbind(sc13a_1, s13aoa_1, sc13j_1, s13joa_1, sc13mt_1, s13moa_1), na.rm = TRUE),
         predSC18 = rowMeans(cbind(sc18a_1, s18aoa_1, sc18j_1, s18joa_1, sc18mt_1, s18moa_1), na.rm = TRUE),
         predSC19 = rowMeans(cbind(sc19a_1, s19aoa_1, sc19j_1, s19joa_1, sc19mt_1, s19moa_1), na.rm = TRUE),
         predbSC2 = rowMeans(cbind(sc2a_2, s2aoa_2, sc2j_2, s2joa_2, sc2mt_2, s2moa_2), na.rm = TRUE),
         predbSC4 = rowMeans(cbind(sc4a_2, s4aoa_2, sc4j_2, s4joa_2, sc4mt_2, s4moa_2), na.rm = TRUE),
         predbSC7 = rowMeans(cbind(sc7a_2, s7aoa_2, sc7j_2, s7joa_2, sc7mt_2, s7moa_2), na.rm = TRUE),
         predbSC11 = rowMeans(cbind(sc11a_2, s11aoa_2, sc11j_2, s11joa_2, sc11mt_2, s11moa_2), na.rm = TRUE),
         predbSC12 = rowMeans(cbind(sc12a_2, s12aoa_2, sc12j_2, s12joa_2, sc12mt_2, s12moa_2), na.rm = TRUE),
         predbSC13 = rowMeans(cbind(sc13a_2, s13aoa_2, sc13j_2, s13joa_2, sc13mt_2, s13moa_2), na.rm = TRUE),
         predbSC18 = rowMeans(cbind(sc18a_2, s18aoa_2, sc18j_2, s18joa_2, sc18mt_2, s18moa_2), na.rm = TRUE),
         predbSC19 = rowMeans(cbind(sc19a_2, s19aoa_2, sc19j_2, s19joa_2, sc19mt_2, s19moa_2), na.rm = TRUE))

cj2 <- cj2 %>%
  mutate(predSC2 = rowMeans(cbind(s2aoa_1, s2joa_1, s2moa_1), na.rm = TRUE),
         predSC4 = rowMeans(cbind(s4aoa_1, s4joa_1, s4moa_1), na.rm = TRUE),
         predSC7 = rowMeans(cbind(s7aoa_1, s7joa_1, s7moa_1), na.rm = TRUE),
         predSC11 = rowMeans(cbind(s11aoa_1, s11joa_1, s11moa_1), na.rm = TRUE),
         predSC12 = rowMeans(cbind(s12aoa_1, s12joa_1, s12moa_1), na.rm = TRUE),
         predSC13 = rowMeans(cbind(s13aoa_1, s13joa_1, s13moa_1), na.rm = TRUE),
         predSC18 = rowMeans(cbind(s18aoa_1, s18joa_1, s18moa_1), na.rm = TRUE),
         predSC19 = rowMeans(cbind(s19aoa_1, s19joa_1, s19moa_1), na.rm = TRUE),
         predbSC2 = rowMeans(cbind(s2aoa_2, s2joa_2, s2moa_2), na.rm = TRUE),
         predbSC4 = rowMeans(cbind(s4aoa_2, s4joa_2, s4moa_2), na.rm = TRUE),
         predbSC7 = rowMeans(cbind(s7aoa_2, s7joa_2, s7moa_2), na.rm = TRUE),
         predbSC11 = rowMeans(cbind(s11aoa_2, s11joa_2, s11moa_2), na.rm = TRUE),
         predbSC12 = rowMeans(cbind(s12aoa_2, s12joa_2, s12moa_2), na.rm = TRUE),
         predbSC13 = rowMeans(cbind(s13aoa_2, s13joa_2, s13moa_2), na.rm = TRUE),
         predbSC18 = rowMeans(cbind(s18aoa_2, s18joa_2, s18moa_2), na.rm = TRUE),
         predbSC19 = rowMeans(cbind(s19aoa_2, s19joa_2, s19moa_2), na.rm = TRUE))

# Record initial predictions in the same way
cj1 <- cj1 %>%
  mutate(initSC2 = rowMeans(cbind(sc2_1, s2ao_1, s2jo_1, s2mo_1), na.rm = TRUE),
         initSC4 = rowMeans(cbind(sc4_1, s4ao_1, s4jo_1, s4mo_1), na.rm = TRUE),
         initSC7 = rowMeans(cbind(sc7_1, s7ao_1, s7jo_1, s7mo_1), na.rm = TRUE),
         initSC11 = rowMeans(cbind(sc11_1, s11ao_1, s11jo_1, s11mo_1), na.rm = TRUE),
         initSC12 = rowMeans(cbind(sc12_1, s12ao_1, s12jo_1, s12mo_1), na.rm = TRUE),
         initSC13 = rowMeans(cbind(sc13_1, s13ao_1, s13jo_1, s13mo_1), na.rm = TRUE),
         initSC18 = rowMeans(cbind(sc18_1, s18ao_1, s18jo_1, s18mo_1), na.rm = TRUE),
         initSC19 = rowMeans(cbind(sc19_1, s19ao_1, s19jo_1, s19mo_1), na.rm = TRUE),
         initbSC2 = rowMeans(cbind(sc2_2, s2ao_2, s2jo_2, s2mo_2), na.rm = TRUE),
         initbSC4 = rowMeans(cbind(sc4_2, s4ao_2, s4jo_2, s4mo_2), na.rm = TRUE),
         initbSC7 = rowMeans(cbind(sc7_2, s7ao_2, s7jo_2, s7mo_2), na.rm = TRUE),
         initbSC11 = rowMeans(cbind(sc11_2, s11ao_2, s11jo_2, s11mo_2), na.rm = TRUE),
         initbSC12 = rowMeans(cbind(sc12_2, s12ao_2, s12jo_2, s12mo_2), na.rm = TRUE),
         initbSC13 = rowMeans(cbind(sc13_2, s13ao_2, s13jo_2, s13mo_2), na.rm = TRUE),
         initbSC18 = rowMeans(cbind(sc18_2, s18ao_2, s18jo_2, s18mo_2), na.rm = TRUE),
         initbSC19 = rowMeans(cbind(sc19_2, s19ao_2, s19jo_2, s19mo_2), na.rm = TRUE))

cj2 <- cj2 %>%
  mutate(initSC2 = rowMeans(cbind(s2ao_1, s2jo_1, s2mo_1), na.rm = TRUE),
         initSC4 = rowMeans(cbind(s4ao_1, s4jo_1, s4mo_1), na.rm = TRUE),
         initSC7 = rowMeans(cbind(s7ao_1, s7jo_1, s7mo_1), na.rm = TRUE),
         initSC11 = rowMeans(cbind(s11ao_1, s11jo_1, s11mo_1), na.rm = TRUE),
         initSC12 = rowMeans(cbind(s12ao_1, s12jo_1, s12mo_1), na.rm = TRUE),
         initSC13 = rowMeans(cbind(s13ao_1, s13jo_1, s13mo_1), na.rm = TRUE),
         initSC18 = rowMeans(cbind(s18ao_1, s18jo_1, s18mo_1), na.rm = TRUE),
         initSC19 = rowMeans(cbind(s19ao_1, s19jo_1, s19mo_1), na.rm = TRUE),
         initbSC2 = rowMeans(cbind(s2ao_2, s2jo_2, s2mo_2), na.rm = TRUE),
         initbSC4 = rowMeans(cbind(s4ao_2, s4jo_2, s4mo_2), na.rm = TRUE),
         initbSC7 = rowMeans(cbind(s7ao_2, s7jo_2, s7mo_2), na.rm = TRUE),
         initbSC11 = rowMeans(cbind(s11ao_2, s11jo_2, s11mo_2), na.rm = TRUE),
         initbSC12 = rowMeans(cbind(s12ao_2, s12jo_2, s12mo_2), na.rm = TRUE),
         initbSC13 = rowMeans(cbind(s13ao_2, s13jo_2, s13mo_2), na.rm = TRUE),
         initbSC18 = rowMeans(cbind(s18ao_2, s18jo_2, s18mo_2), na.rm = TRUE),
         initbSC19 = rowMeans(cbind(s19ao_2, s19jo_2, s19mo_2), na.rm = TRUE))

# Break up the data by scenario to stack the data
sc2 <- cj1 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female,
         partisanship, needcognition, needjudge, edvalue, ends_with("SC2"), source) %>% 
  mutate(scenario = 2)
names(sc2) <- c("ResponseId", "tia", "age", "encouragement", "female",
                "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                "initb","source","scenario")
sc4 <- cj1 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC4"), source) %>%
  mutate(scenario = 4)
names(sc4) <- c("ResponseId", "tia", "age", "encouragement", "female",
                "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                "initb","source","scenario")
sc7 <- cj1 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC7"), source) %>%
  mutate(scenario = 7)
names(sc7) <- c("ResponseId", "tia", "age", "encouragement", "female",
                "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                "initb","source","scenario")
sc11 <- cj1 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC11"), source) %>%
  mutate(scenario = 11)
names(sc11) <- c("ResponseId", "tia", "age", "encouragement", "female",
                 "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                "initb","source","scenario")
sc12 <- cj1 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC12"), source) %>%
  mutate(scenario = 12)
names(sc12) <- c("ResponseId", "tia", "age", "encouragement", "female",
                 "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                "initb","source","scenario")
sc13 <- cj1 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC13"), source) %>%
  mutate(scenario = 13)
names(sc13) <- c("ResponseId", "tia", "age", "encouragement", "female",
                 "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                "initb","source","scenario")
sc18 <- cj1 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC18"), source) %>%
  mutate(scenario = 18)
names(sc18) <- c("ResponseId", "tia", "age", "encouragement", "female",
                 "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                "initb","source","scenario")
sc19 <- cj1 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC19"), source) %>%
  mutate(scenario = 19)
names(sc19) <- c("ResponseId", "tia", "age", "encouragement", "female",
                 "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                "initb","source","scenario")

sc2a <- cj2 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC2"), source) %>% 
  mutate(scenario = 2)
names(sc2a) <- c("ResponseId", "tia", "age", "encouragement", "female",
                 "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                 "initb","source","scenario")
sc4a <- cj2 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC4"), source) %>% 
  mutate(scenario = 4)
names(sc4a) <- c("ResponseId", "tia", "age", "encouragement", "female",
                 "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                 "initb","source","scenario")
sc7a <- cj2 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC7"), source) %>% 
  mutate(scenario = 7)
names(sc7a) <- c("ResponseId", "tia", "age", "encouragement", "female",
                 "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                 "initb","source","scenario")
sc11a <- cj2 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC11"), source) %>% 
  mutate(scenario = 11)
names(sc11a) <- c("ResponseId", "tia", "age", "encouragement", "female",
                  "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                 "initb","source","scenario")
sc12a <- cj2 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC12"), source) %>% 
  mutate(scenario = 12)
names(sc12a) <- c("ResponseId", "tia", "age", "encouragement", "female",
                  "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                 "initb","source","scenario")
sc13a <- cj2 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC13"), source) %>% 
  mutate(scenario = 13)
names(sc13a) <- c("ResponseId", "tia", "age", "encouragement", "female",
                  "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                 "initb","source","scenario")
sc18a <- cj2 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC18"), source) %>% 
  mutate(scenario = 18)
names(sc18a) <- c("ResponseId", "tia", "age", "encouragement", "female",
                  "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                 "initb","source","scenario")
sc19a <- cj2 %>% 
  select(ResponseId, trustAutomation, age, encouragement, female, 
         partisanship, needcognition, needjudge, edvalue, ends_with("SC19"), source) %>% 
  mutate(scenario = 19)
names(sc19a) <- c("ResponseId", "tia", "age", "encouragement", "female",
                 "partisanship","needcog", "needjudge", "ed","advice", "adviceb", "adType", "pred", "predb", "init", 
                 "initb","source","scenario")


cj1long <- rbind(sc2, sc4, sc7, sc11, sc12, sc13, sc18, sc19,
                 sc2a, sc4a, sc7a, sc11a, sc12a, sc13a, sc18a, sc19a)

cj1long <- cj1long %>%
  mutate(distance = (abs(pred - advice) + abs(predb - adviceb))/2,
         adviceWt = ((abs(pred - init)/abs(advice - init)) + (abs(predb - initb)/abs(adviceb - initb)))/2,
         brier = (((0.01*advice) - (0.01*pred))^2 + ((0.01*adviceb) - (0.01*predb))^2)/2,
         algorithm = ifelse(adType == "algorithm", 1, 0),
         judge = ifelse(adType == "judge", 1, 0))


#### Table of summary statistics
# Table A2 in SI
experimentSum <- cj1long %>%
  select(distance, adviceWt, brier, algorithm, judge,
         tia, age, ed, encouragement, female, partisanship, 
         needcog, needjudge) %>%
  mutate(adviceWt = ifelse(adviceWt > 1, NA, adviceWt))

experimentSum <- data.frame(experimentSum)

sink(here("cjSummaries.tex"))
stargazer(experimentSum, style="ajps", 
          title="Summary Statistics", 
          covariate.labels=c("Distance to Advice", "Weight of Advice", "Squared Distance",
                             "Algorithm Treatment","Judge Treatment",
                             "Trust in Automation", "Age",
                             "Education", "Encouragement", "Female",
                             "Partisanship", "Need for Cognition", "Need to Evaluate")
)
sink()



#### Graphs of baseline weight of advice
# Figure A3 in SI
cj1long %>%
  filter(adviceWt <= 1) %>%
  mutate(adType = case_when(adType == "algorithm" ~ "Algorithm",
                            adType == "judge" ~ "Judge",
                            adType == "mt" ~ "Previous Respondents")) %>%
  ggplot() + geom_histogram(aes(x = adviceWt, fill = adType), binwidth = .1) + 
  facet_wrap(~adType, ncol = 2) + theme_bw() +
  labs(x = "Advice Weight", y = "Count", fill = "Source")
# Summary statistics for baselines as reported in text
cj1long %>%
  filter(adviceWt <= 1) %>%
  summarize(avg = mean(adviceWt, na.rm = T),
            med = median(adviceWt, na.rm = T)) 
cj1long %>%
  filter(adviceWt <= 1 & source == "Lucid") %>%
  summarize(avg = mean(adviceWt, na.rm = T),
            med = median(adviceWt, na.rm = T)) 
cj1long %>%
  filter(adviceWt <= 1 & source == "MTurk") %>%
  summarize(avg = mean(adviceWt, na.rm = T),
            med = median(adviceWt, na.rm = T)) 
# Figure A4 in SI
acc1 <- cj1long %>%
  filter(adviceWt <= 1) %>%
  mutate(acceptAdvice = ifelse(adviceWt > 0.5, 1, 0)) %>%
  group_by(adType, source) %>%
  summarise(respondents = n(),
            sumAccept = sum(acceptAdvice)) %>%
  ungroup() %>%
  mutate(pctAccept = sumAccept/respondents, 
         adType = case_when(adType == "algorithm" ~ "Algorithm",
                            adType == "judge" ~ "Judge",
                            adType == "mt" ~ "Previous Respondents")) %>%
  ggplot() + geom_bar(aes(x = adType, y = pctAccept, fill = source), stat = "identity", position = "dodge") +
  scale_y_continuous(breaks = seq(0, 0.3, 0.1)) +
  xlab("Source Treatment") + ylab("Percent Accepting Advice (Weight > 0.5)") +
  theme_bw()
acc1



#### Balance models for SI
# Table A6 in SI
bal1 <- lmer(algorithm ~ age + ed + female + partisanship +
               (1|ResponseId), 
             data = cj1long[cj1long$adviceWt <= 1,],
             control = lmerControl(optimizer = "Nelder_Mead"))
summary(bal1)

bal2 <- lmer(judge ~ age + ed + female + partisanship +
               (1|ResponseId), 
             data = cj1long[cj1long$adviceWt <= 1,],
             control = lmerControl(optimizer = "Nelder_Mead"))
summary(bal2)

# Balance table
# Function for adding rows to a data frame (will be used later)
insertrow <- function(existingDF, newrow, r) {
  existingDF[seq(r+1,nrow(existingDF)+1),] <- existingDF[seq(r,nrow(existingDF)),]
  existingDF[r,] <- newrow
  existingDF
}
# Create a baseline table
Tables <- stargazer(bal1, bal2, style="ajps", 
                    title="Criminal Justice Balance Table", 
                    dep.var.labels.include = FALSE, 
                    covariate.labels=c("Age","Education","Female", "Republican Partisanship")
)
# Convert Tables to data frame and coerce the one vector (also called Tables) to character
Tables <- as.data.frame(Tables)
Tables$Tables <- as.character(Tables$Tables)
# Find where you want to put in the random effect. In our case, this is right after the last fixed effect. Line: 23.
r <- 21
# Create some standard label lines
randomeffect <- "{\\bf Random Effect} & \\\\"
hline <- "\\hline"
newline <- "\\\\"
# Now insert those label lines where they are needed
Tables <- insertrow(Tables, hline, r)
Tables <- insertrow(Tables,randomeffect,r+1)
Tables <- insertrow(Tables,hline,r+2)
# Get number of unique values in each grouping
num.participants <- sapply(ranef(bal1),nrow)[1]
# Get standard deviation of the random effect
stddev.model1.participants <- attributes(VarCorr(bal1)$ResponseId)$stddev
stddev.model2.participants <- attributes(VarCorr(bal2)$ResponseId)$stddev
# Create a LaTex character row for the random effect
number.of.participants <- paste("\\# of Participants & ", num.participants, "&", num.participants, "\\\\")
stddev.participants <- paste("Participant Standard Deviation & ", round(stddev.model1.participants, 3), "&", round(stddev.model2.participants, 3), "\\\\")
# Add these lines to the table
Tables <- insertrow(Tables,number.of.participants,r+3)
Tables <- insertrow(Tables,stddev.participants,r+4)
Tables <- insertrow(Tables,newline,r+5)
# Write the table to a file that can be inserted into the document
write.table(Tables,file=here("CJBalance.tex"),
            sep="",row.names= FALSE,na="", quote = FALSE, col.names = FALSE)


#### Main models
# Figure 2 and Table A9
model1b <- lmer(distance ~ algorithm + judge + (1|ResponseId) + (1|scenario), 
                data = cj1long, control = lmerControl(optimizer = "Nelder_Mead"))
summary(model1b)
model2b <- lmer(adviceWt ~ algorithm + judge + (1|ResponseId) + (1|scenario), 
                data = cj1long[cj1long$adviceWt <= 1,],
                control = lmerControl(optimizer = "Nelder_Mead"))
summary(model2b)
# Plots for Figure 2
# Model 1
mod1sum <- summary(model1b)
mod1sim <- sim(model1b, n.sims = 1000)
mod1simcoef <- as_tibble(data.frame(coef(mod1sim)$fixef))
mod1simcoef$simulation <- rownames(mod1simcoef)
mod1simcoeflong <- gather(mod1simcoef, condition, coefficient, X.Intercept.:judge)
mod1simcoeflong$condition[mod1simcoeflong$condition == "algorithm"] <- "Computer Algorithm"
mod1simcoeflong$condition[mod1simcoeflong$condition == "judge"] <- "Judge"

m1b95 <- mod1simcoeflong %>%
  filter(condition == "Computer Algorithm" | condition == "Judge") %>%
  group_by(condition) %>%
  summarize(mcond = mean(coefficient, na.rm = T),
            sdcond = sd(coefficient, na.rm = T)) %>%
  mutate(hi = mcond + 1.96*sdcond,
         lo = mcond - 1.96*sdcond) %>% 
  ggplot() +
  geom_point(aes(x = condition, y = mcond), size = 3, color = "red") +
  geom_errorbar(aes(x = condition, ymin = lo, ymax = hi), color = "red") +
  geom_hline(aes(yintercept = 0)) +
  labs(x = "Condition", y = "ATE", title = "Distance to Advice") +
  coord_flip() + theme_bw()
m1b95

# Model 2
mod2sum <- summary(model2b)
mod2sim <- sim(model2b, n.sims = 1000)
mod2simcoef <- as_tibble(data.frame(coef(mod2sim)$fixef))
mod2simcoef$simulation <- rownames(mod2simcoef)
mod2simcoeflong <- gather(mod2simcoef, condition, coefficient, X.Intercept.:judge)
mod2simcoeflong$condition[mod2simcoeflong$condition == "algorithm"] <- "Computer Algorithm"
mod2simcoeflong$condition[mod2simcoeflong$condition == "judge"] <- "Judge"

m2b95 <- mod2simcoeflong %>%
  filter(condition == "Computer Algorithm" | condition == "Judge") %>%
  group_by(condition) %>%
  summarize(mcond = mean(coefficient, na.rm = T),
            sdcond = sd(coefficient, na.rm = T)) %>%
  mutate(hi = mcond + 1.96*sdcond,
         lo = mcond - 1.96*sdcond) %>% 
  ggplot() +
  geom_point(aes(x = condition, y = mcond), size = 3, color = "blue") +
  geom_errorbar(aes(x = condition, ymin = lo, ymax = hi), color = "blue") +
  geom_hline(aes(yintercept = 0)) +
  labs(x = "Condition", y = "ATE", title = "Weight of Advice") +
  coord_flip() + theme_bw()
m2b95

g1 <- grid.arrange(m1b95, m2b95, ncol = 1)

# Create Table A9
# Function for adding rows to a data frame (will be used later)
insertrow <- function(existingDF, newrow, r) {
  existingDF[seq(r + 1,nrow(existingDF) + 1),] <- existingDF[seq(r,nrow(existingDF)),]
  existingDF[r,] <- newrow
  existingDF
}
# Create a baseline table
Tables <- stargazer(model1b, model2b, style = "ajps", 
                    title = "Models for Weight of Advice", 
                    dep.var.labels.include = FALSE, 
                    covariate.labels = c("Algorithm", "Judge")
)
# Convert Tables to data frame and coerce the one vector (also called Tables) to character
Tables <- as.data.frame(Tables)
Tables$Tables <- as.character(Tables$Tables)
# Find where you want to put in the random effect. In our case, this is right after the last fixed effect. Line: 23.
r <- 15
# Create some standard label lines
randomeffect <- "{\\bf Random Effect} & \\\\"
hline <- "\\hline"
newline <- "\\\\"
# Now insert those label lines where they are needed
Tables <- insertrow(Tables, hline, r)
Tables <- insertrow(Tables,randomeffect,r + 1)
Tables <- insertrow(Tables,hline,r + 2)
# Get number of unique values in each grouping
num.participants1 <- sapply(ranef(model1b),nrow)[1]
num.scenarios1 <- sapply(ranef(model1b),nrow)[2]
num.participants2 <- sapply(ranef(model2b),nrow)[1]
num.scenarios2 <- sapply(ranef(model2b),nrow)[2]
# Get standard deviation of the random effect
stddev.model1.participants <- attributes(VarCorr(model1b)$ResponseId)$stddev
stddev.model1.scenarios <- attributes(VarCorr(model1b)$scenario)$stddev
stddev.model2.participants <- attributes(VarCorr(model2b)$ResponseId)$stddev
stddev.model2.scenarios <- attributes(VarCorr(model2b)$scenario)$stddev
# Create a LaTex character row for the random effect
number.of.participants <- paste("\\# of Participants & ", num.participants1, " & ", num.participants2, "\\\\")
stddev.participants <- paste("Participant Standard Deviation & ", round(stddev.model1.participants, 3), " & ", round(stddev.model2.participants, 3), "\\\\")
number.of.scenarios <- paste("\\# of Scenarios & ", num.scenarios1, " & ", num.scenarios2, "\\\\")
stddev.scenarios <- paste("Scenario Standard Deviation & ", round(stddev.model1.scenarios, 3), " & ", round(stddev.model2.scenarios, 3), "\\\\")
# Add these lines to the table
Tables <- insertrow(Tables,number.of.participants,r + 3)
Tables <- insertrow(Tables,stddev.participants,r + 4)
Tables <- insertrow(Tables,newline,r + 5)
Tables <- insertrow(Tables,number.of.scenarios,r + 6)
Tables <- insertrow(Tables,stddev.scenarios,r + 7)
Tables <- insertrow(Tables, hline, r + 8)
# Write the table to a file that can be inserted into the document
write.table(Tables,file = here("REtable(NC).tex"),
            sep = "", row.names = FALSE, na = "", 
            quote = FALSE, col.names = FALSE)


#### Models with alternative baseline
#Table A13 in SI
cjlongb <- cj1long %>%
  mutate(avghuman = ifelse(algorithm == 0 & judge == 0, 1, 0))

model1c <- lmer(distance ~ algorithm + avghuman + (1|ResponseId) + (1|scenario), 
                data = cjlongb, control = lmerControl(optimizer = "Nelder_Mead"))
summary(model1c)
model2c <- lmer(adviceWt ~ algorithm + avghuman + (1|ResponseId) + (1|scenario), 
                data = cjlongb[cjlongb$adviceWt <= 1,],
                control = lmerControl(optimizer = "Nelder_Mead"))
summary(model2c)
model3c <- lmer(brier ~ algorithm + avghuman + (1|ResponseId) + (1|scenario), 
                data = cjlongb, control = lmerControl(optimizer = "Nelder_Mead"))
summary(model3c)
# Table generation
# Function for adding rows to a data frame (will be used later)
insertrow <- function(existingDF, newrow, r) {
  existingDF[seq(r + 1,nrow(existingDF) + 1),] <- existingDF[seq(r,nrow(existingDF)),]
  existingDF[r,] <- newrow
  existingDF
}
# Create a baseline table
Tables <- stargazer(model1c, model2c, model3c, style = "ajps", 
                    title = "Models for Weight of Advice", 
                    dep.var.labels.include = FALSE, 
                    covariate.labels = c("Algorithm", "Average Human")
)
# Convert Tables to data frame and coerce the one vector (also called Tables) to character
Tables <- as.data.frame(Tables)
Tables$Tables <- as.character(Tables$Tables)
# Find where you want to put in the random effect. In our case, this is right after the last fixed effect. Line: 23.
r <- 15
# Create some standard label lines
randomeffect <- "{\\bf Random Effect} & \\\\"
hline <- "\\hline"
newline <- "\\\\"
# Now insert those label lines where they are needed
Tables <- insertrow(Tables, hline, r)
Tables <- insertrow(Tables,randomeffect,r + 1)
Tables <- insertrow(Tables,hline,r + 2)
# Get number of unique values in each grouping
num.participants1 <- sapply(ranef(model1c),nrow)[1]
num.scenarios1 <- sapply(ranef(model1c),nrow)[2]
num.participants2 <- sapply(ranef(model2c),nrow)[1]
num.scenarios2 <- sapply(ranef(model2c),nrow)[2]
num.participants3 <- sapply(ranef(model3c),nrow)[1]
num.scenarios3 <- sapply(ranef(model3c),nrow)[2]
# Get standard deviation of the random effect
stddev.model1.participants <- attributes(VarCorr(model1c)$ResponseId)$stddev
stddev.model1.scenarios <- attributes(VarCorr(model1c)$scenario)$stddev
stddev.model2.participants <- attributes(VarCorr(model2c)$ResponseId)$stddev
stddev.model2.scenarios <- attributes(VarCorr(model2c)$scenario)$stddev
stddev.model3.participants <- attributes(VarCorr(model3c)$ResponseId)$stddev
stddev.model3.scenarios <- attributes(VarCorr(model3c)$scenario)$stddev
# Create a LaTex character row for the random effect
number.of.participants <- paste("\\# of Participants & ", num.participants1, " & ", num.participants2, " & ", num.participants3, "\\\\")
stddev.participants <- paste("Participant Standard Deviation & ", round(stddev.model1.participants, 3), " & ", round(stddev.model2.participants, 3), " & ", round(stddev.model3.participants, 3), "\\\\")
number.of.scenarios <- paste("\\# of Scenarios & ", num.scenarios1, " & ", num.scenarios2, " & ", num.scenarios3, "\\\\")
stddev.scenarios <- paste("Scenario Standard Deviation & ", round(stddev.model1.scenarios, 3), " & ", round(stddev.model2.scenarios, 3), " & ", round(stddev.model3.scenarios, 3), "\\\\")
# Add these lines to the table
Tables <- insertrow(Tables,number.of.participants,r + 3)
Tables <- insertrow(Tables,stddev.participants,r + 4)
Tables <- insertrow(Tables,newline,r + 5)
Tables <- insertrow(Tables,number.of.scenarios,r + 6)
Tables <- insertrow(Tables,stddev.scenarios,r + 7)
Tables <- insertrow(Tables, hline, r + 8)
# Write the table to a file that can be inserted into the document
write.table(Tables,file = here("REtable(AltB).tex"),
            sep = "", row.names = FALSE, na = "", 
            quote = FALSE, col.names = FALSE)



#### Models with controls
# Models for tia with controls
model1 <- lmer(distance ~ algorithm + judge + tia + age + ed + encouragement + 
                 female + partisanship + needcog + needjudge +
                 (1|ResponseId) + (1|scenario), 
               data = cj1long[cj1long$adviceWt <= 1,],
               control = lmerControl(optimizer = "Nelder_Mead"))
summary(model1)
model2 <- lmer(adviceWt ~ algorithm + judge + tia + age + ed + encouragement + 
                 female + partisanship + needcog + needjudge +
                 (1|ResponseId) + (1|scenario), 
               data = cj1long[cj1long$adviceWt <= 1,],
               control = lmerControl(optimizer = "Nelder_Mead"))
summary(model2)
# Generate Table A17
# Function for adding rows to a data frame (will be used later)
insertrow <- function(existingDF, newrow, r) {
  existingDF[seq(r+1,nrow(existingDF)+1),] <- existingDF[seq(r,nrow(existingDF)),]
  existingDF[r,] <- newrow
  existingDF
}
# Create a baseline table
Tables <- stargazer(model1, model2, style="ajps", 
                    title="Models for Distance and Weight of Advice", 
                    dep.var.labels.include = FALSE, 
                    covariate.labels=c( "Algorithm", "Judge", "Trust in Automation", 
                                        "Age","Education","Encouragement", "Female", "Republican Partisanship", "Need for Cognition",
                                        "Need to Evaluate")
)
# Convert Tables to data frame and coerce the one vector (also called Tables) to character
Tables <- as.data.frame(Tables)
Tables$Tables <- as.character(Tables$Tables)
# Find where you want to put in the random effect. In our case, this is right after the last fixed effect. Line: 23.
r <- 31
# Create some standard label lines
randomeffect <- "{\\bf Random Effect} & \\\\"
hline <- "\\hline"
newline <- "\\\\"
# Now insert those label lines where they are needed
Tables <- insertrow(Tables, hline, r)
Tables <- insertrow(Tables,randomeffect,r+1)
Tables <- insertrow(Tables,hline,r+2)
# Get number of unique values in each grouping
num.participants1 <- sapply(ranef(model1),nrow)[1]
num.scenarios1 <- sapply(ranef(model1),nrow)[2]
num.participants2 <- sapply(ranef(model2),nrow)[1]
num.scenarios2 <- sapply(ranef(model2),nrow)[2]
# Get standard deviation of the random effect
stddev.model1.participants <- attributes(VarCorr(model1)$ResponseId)$stddev
stddev.model1.scenarios <- attributes(VarCorr(model1)$scenario)$stddev
stddev.model2.participants <- attributes(VarCorr(model2)$ResponseId)$stddev
stddev.model2.scenarios <- attributes(VarCorr(model2)$scenario)$stddev
# Create a LaTex character row for the random effect
number.of.participants <- paste("\\# of Participants & ", num.participants1, " & ", num.participants2, "\\\\")
stddev.participants <- paste("Participant Standard Deviation & ", round(stddev.model1.participants, 3), " & ", round(stddev.model2.participants, 3), "\\\\")
number.of.scenarios <- paste("\\# of Scenarios & ", num.scenarios1, " & ", num.scenarios2, "\\\\")
stddev.scenarios <- paste("Scenario Standard Deviation & ", round(stddev.model1.scenarios, 3), " & ", round(stddev.model2.scenarios, 3), "\\\\")
# Add these lines to the table
Tables <- insertrow(Tables,number.of.participants,r+3)
Tables <- insertrow(Tables,stddev.participants,r+4)
Tables <- insertrow(Tables,newline,r+5)
Tables <- insertrow(Tables,number.of.scenarios,r+6)
Tables <- insertrow(Tables,stddev.scenarios,r+7)
Tables <- insertrow(Tables, hline, r+8)
# Write the table to a file that can be inserted into the document
write.table(Tables,file=here("REtable(wcov)cj.tex"),
            sep="",row.names= FALSE,na="", quote = FALSE, col.names = FALSE)



#### Heterogeneous treatment effects
# Right side of Figure A22 and Table A15 in SI
model5 <- lmer(adviceWt ~ algorithm + judge + I(algorithm*tia) + tia + age + ed + encouragement + 
                 female + partisanship + needcog + needjudge +
                 (1|ResponseId) + (1|scenario), 
               data = cj1long[cj1long$adviceWt <= 1,],
               control = lmerControl(optimizer = "Nelder_Mead"))
summary(model5)
model8 <- lmer(adviceWt ~ algorithm + judge + I(algorithm*age) + tia + age + ed + encouragement + 
                 female + partisanship + needcog + needjudge +
                 (1|ResponseId) + (1|scenario), 
               data = cj1long[cj1long$adviceWt <= 1,],
               control = lmerControl(optimizer = "Nelder_Mead"))
summary(model8)
model11 <- lmer(adviceWt ~ algorithm + judge + I(algorithm*ed) + tia + age + ed + encouragement + 
                  female + partisanship + needcog + needjudge +
                  (1|ResponseId) + (1|scenario), 
                data = cj1long[cj1long$adviceWt <= 1,],
                control = lmerControl(optimizer = "Nelder_Mead"))
summary(model11)
model14 <- lmer(adviceWt ~ algorithm + judge + I(algorithm*female) + tia + age + ed + encouragement + 
                  female + partisanship + needcog + needjudge +
                  (1|ResponseId) + (1|scenario), 
                data = cj1long[cj1long$adviceWt <= 1,])
summary(model14)

# Model 5 Interaction Plot
mod5sum <- summary(model5)
tiaX <- seq(1,7)
maineffect5 <- fixef(model5)[2] + fixef(model5)[4] * tiaX
maineffect5 <- as_tibble(cbind(tiaX,maineffect5))
mod5sim <- sim(model5, n.sims = 1000)
mod5simcoef <- as.tibble(data.frame(coef(mod5sim)$fixef))
mod5simcoef$simulation <- rownames(mod5simcoef)
mod5simcoeflong <- mod5simcoef %>%
  mutate(simeffect1 = algorithm + I.algorithm...tia. * 1,
         simeffect2 = algorithm + I.algorithm...tia. * 2,
         simeffect3 = algorithm + I.algorithm...tia. * 3,
         simeffect4 = algorithm + I.algorithm...tia. * 4,
         simeffect5 = algorithm + I.algorithm...tia. * 5,
         simeffect6 = algorithm + I.algorithm...tia. * 6,
         simeffect7 = algorithm + I.algorithm...tia. * 7) %>%
  dplyr::select(simulation:simeffect7) %>%
  gather(tialevel, estimate, simeffect1:simeffect7) %>%
  mutate(tia = as.numeric(substring(tialevel, nchar(tialevel), nchar(tialevel))))

p4 = ggplot() + geom_line(data = mod5simcoeflong,aes(x = tia,y = estimate, group = simulation), color = "blue", alpha = 0.1) +
  geom_line(data = maineffect5, aes(x = tiaX, y = maineffect5), size = 1.5) +
  xlab("Trust in Automation") + ylab("Estimated Impact of Algorithm") + ggtitle("Weight of Advice") +
  geom_hline(aes(yintercept = 0)) + theme_bw()
p4

# Model 8 Interaction Plot
mod8sum <- summary(model8)
tiaX <- seq(19,74, by = 5)
maineffect8 <- fixef(model8)[2] + fixef(model8)[4] * tiaX
maineffect8 <- as_tibble(cbind(tiaX,maineffect8))
mod8sim <- sim(model8, n.sims = 1000)
mod8simcoef <- as.tibble(data.frame(coef(mod8sim)$fixef))
mod8simcoef$simulation <- rownames(mod8simcoef)
mod8simcoeflong <- mod8simcoef %>%
  mutate(simeffect19 = algorithm + I.algorithm...age. * 19,
         simeffect24 = algorithm + I.algorithm...age. * 24,
         simeffect29 = algorithm + I.algorithm...age. * 29,
         simeffect34 = algorithm + I.algorithm...age. * 34,
         simeffect39 = algorithm + I.algorithm...age. * 39,
         simeffect44 = algorithm + I.algorithm...age. * 44,
         simeffect49 = algorithm + I.algorithm...age. * 49,
         simeffect54 = algorithm + I.algorithm...age. * 54,
         simeffect59 = algorithm + I.algorithm...age. * 59,
         simeffect64 = algorithm + I.algorithm...age. * 64,
         simeffect69 = algorithm + I.algorithm...age. * 69,
         simeffect74 = algorithm + I.algorithm...age. * 74) %>%
  dplyr::select(simulation:simeffect74) %>%
  gather(tialevel, estimate, simeffect19:simeffect74) %>%
  mutate(tia = as.numeric(substring(tialevel, nchar(tialevel) - 1, nchar(tialevel))))

p6 = ggplot() + geom_line(data = mod8simcoeflong,aes(x = tia,y = estimate, group = simulation), color = "blue", alpha = 0.1) +
  geom_line(data = maineffect8, aes(x = tiaX, y = maineffect8), size = 1.5) +
  xlab("Age") + ylab("Estimated Impact of Algorithm") + ggtitle("Weight of Advice") +
  geom_hline(aes(yintercept = 0)) + theme_bw()
p6

# Interaction graph for model 11
mod11sum <- summary(model11)
tiaX <- seq(2,4)
maineffect11 <- fixef(model11)[2] + fixef(model11)[4] * tiaX
maineffect11 <- as_tibble(cbind(tiaX,maineffect11))
mod11sim <- sim(model11, n.sims = 1000)
mod11simcoef <- as.tibble(data.frame(coef(mod11sim)$fixef))
mod11simcoef$simulation <- rownames(mod11simcoef)
mod11simcoeflong <- mod11simcoef %>%
  mutate(simeffect2 = algorithm + I.algorithm...ed. * 2,
         simeffect3 = algorithm + I.algorithm...ed. * 3,
         simeffect4 = algorithm + I.algorithm...ed. * 4) %>%
  dplyr::select(simulation:simeffect4) %>%
  gather(tialevel, estimate, simeffect2:simeffect4) %>%
  mutate(tia = as.numeric(substring(tialevel, nchar(tialevel), nchar(tialevel))))

p8 = ggplot() + geom_line(data = mod11simcoeflong,aes(x = tia,y = estimate, group = simulation), color = "blue", alpha = 0.1) +
  geom_line(data = maineffect11, aes(x = tiaX, y = maineffect11), size = 1.5) +
  xlab("Education") + ylab("Estimated Impact of Algorithm") + ggtitle("Weight of Advice") +
  geom_hline(aes(yintercept = 0)) + theme_bw()
p8

# Interaction graph for model 14
mod14sum <- summary(model14)
tiaX <- c(0,1)
maineffect14 <- fixef(model14)[2] + fixef(model14)[4] * tiaX
maineffect14 <- as_tibble(cbind(tiaX,maineffect14))
mod14sim <- sim(model14, n.sims = 1000)
mod14simcoef <- as.tibble(data.frame(coef(mod14sim)$fixef))
mod14simcoef$simulation <- rownames(mod14simcoef)
mod14simcoeflong <- mod14simcoef %>%
  mutate(simeffect0 = algorithm + I.algorithm...female. * 0,
         simeffect1 = algorithm + I.algorithm...female. * 1) %>%
  dplyr::select(simulation:simeffect1) %>%
  gather(tialevel, estimate, simeffect0:simeffect1) %>%
  mutate(tia = as.numeric(substring(tialevel, nchar(tialevel), nchar(tialevel))))

p10 = ggplot() + geom_line(data = mod14simcoeflong,aes(x = tia,y = estimate, group = simulation), color = "blue", alpha = 0.1) +
  geom_line(data = maineffect14, aes(x = tiaX, y = maineffect14), size = 1.5) +
  xlab("Female") + ylab("Estimated Impact of Algorithm") + ggtitle("Weight of Advice") +
  geom_hline(aes(yintercept = 0)) + theme_bw()
p10

# combine graphs for direct model coefficients
g2b <- grid.arrange(p4,p6,p8,p10,ncol = 1)
ggsave(here(".." ,"..", "Manuscript", "Figures", "cjWhoTrust(weight).pdf"),
       plot = g2b, device = "pdf", width = 4.25, height = 11)



#### Separate models for each wave of criminal justice experiments
# Table A21 in SI
model1l <- lmer(distance ~ algorithm + judge + tia + age + ed + encouragement + 
                  female + partisanship + needcog + needjudge +
                  (1|ResponseId) + (1|scenario), 
                data = cj1long[cj1long$adviceWt <= 1 & cj1long$source == "Lucid",],
                control = lmerControl(optimizer = "Nelder_Mead"))
summary(model1l)
model2l <- lmer(adviceWt ~ algorithm + judge + tia + age + ed + encouragement + 
                  female + partisanship + needcog + needjudge +
                  (1|ResponseId) + (1|scenario), 
                data = cj1long[cj1long$adviceWt <= 1 & cj1long$source == "Lucid",],
                control = lmerControl(optimizer = "Nelder_Mead"))
summary(model2l)
model3l <- lmer(brier ~ algorithm + judge + tia + age + ed + encouragement + 
                  female + partisanship + needcog + needjudge +
                  (1|ResponseId) + (1|scenario), 
                data = cj1long[cj1long$adviceWt <= 1 & cj1long$source == "Lucid",],
                control = lmerControl(optimizer = "Nelder_Mead"))
summary(model3l)
model1m <- lmer(distance ~ algorithm + judge + tia + age + ed + encouragement + 
                  female + partisanship + needcog + needjudge +
                  (1|ResponseId) + (1|scenario), 
                data = cj1long[cj1long$adviceWt <= 1 & cj1long$source == "MTurk",],
                control = lmerControl(optimizer = "Nelder_Mead"))
summary(model1m)
model2m <- lmer(adviceWt ~ algorithm + judge + tia + age + ed + encouragement + 
                  female + partisanship + needcog + needjudge +
                  (1|ResponseId) + (1|scenario), 
                data = cj1long[cj1long$adviceWt <= 1 & cj1long$source == "MTurk",],
                control = lmerControl(optimizer = "Nelder_Mead"))
summary(model2m)
model3m <- lmer(brier ~ algorithm + judge + tia + age + ed + encouragement + 
                  female + partisanship + needcog + needjudge +
                  (1|ResponseId) + (1|scenario), 
                data = cj1long[cj1long$adviceWt <= 1 & cj1long$source == "MTurk",],
                control = lmerControl(optimizer = "Nelder_Mead"))
summary(model3m)

# Table for models separated by source
# Create a baseline table
Tables <- stargazer(model1l, model2l, model3l, model1m, model2m, model3m, style="ajps", 
                    title="Models of Trust in Automation", 
                    dep.var.labels.include = FALSE#, 
)
# Convert Tables to data frame and coerce the one vector (also called Tables) to character
Tables <- as.data.frame(Tables)
Tables$Tables <- as.character(Tables$Tables)
# Find where you want to put in the random effect. In our case, this is right after the last fixed effect. Line: 23.
r <- 31
# Create some standard label lines
randomeffect <- "{\\bf Random Effect} & & \\\\"
hline <- "\\hline"
newline <- "\\\\"
# Now insert those label lines where they are needed
Tables <- insertrow(Tables, hline, r)
Tables <- insertrow(Tables,randomeffect,r+1)
Tables <- insertrow(Tables,hline,r+2)
# Get number of unique values in each grouping
num.participants <- sapply(ranef(model1l),nrow)[1]
num.scenarios <- sapply(ranef(model1l),nrow)[2]
# Get standard deviation of the random effect
stddev.model1.participants <- attributes(VarCorr(model1l)$ResponseId)$stddev
stddev.model1.scenarios <- attributes(VarCorr(model1l)$scenario)$stddev
stddev.model2.participants <- attributes(VarCorr(model2l)$ResponseId)$stddev
stddev.model2.scenarios <- attributes(VarCorr(model2l)$scenario)$stddev
stddev.model3.participants <- attributes(VarCorr(model3l)$ResponseId)$stddev
stddev.model3.scenarios <- attributes(VarCorr(model3l)$scenario)$stddev
stddev.model4.participants <- attributes(VarCorr(model1m)$ResponseId)$stddev
stddev.model4.scenarios <- attributes(VarCorr(model1m)$scenario)$stddev
stddev.model5.participants <- attributes(VarCorr(model2m)$ResponseId)$stddev
stddev.model5.scenarios <- attributes(VarCorr(model2m)$scenario)$stddev
stddev.model6.participants <- attributes(VarCorr(model3m)$ResponseId)$stddev
stddev.model6.scenarios <- attributes(VarCorr(model3m)$scenario)$stddev
# Create a LaTex character row for the random effect
number.of.participants <- paste("\\# of Participants & ", num.participants, "&", num.participants, "&", num.participants, "&", num.participants, "&", num.participants, "&", num.participants, "\\\\")
stddev.participants <- paste("Participant Standard Deviation & ", round(stddev.model1.participants, 3), "&", round(stddev.model2.participants, 3), "&", round(stddev.model3.participants, 3), "&", round(stddev.model4.participants, 3), "&", round(stddev.model5.participants, 3), "&", round(stddev.model6.participants, 3), "\\\\")
number.of.scenarios <- paste("\\# of Scenarios & ", num.scenarios, "&", num.scenarios, "&", num.scenarios, "&", num.scenarios, "&", num.scenarios, "&", num.scenarios, "\\\\")
stddev.scenarios <- paste("Scenario Standard Deviation & ", round(stddev.model1.scenarios, 3), "&", round(stddev.model2.scenarios, 3), "&", round(stddev.model3.scenarios, 3), "&", round(stddev.model4.scenarios, 3), "&", round(stddev.model5.scenarios, 3), "&", round(stddev.model6.scenarios, 3), "\\\\")
# Add these lines to the table
Tables <- insertrow(Tables,number.of.participants,r+3)
Tables <- insertrow(Tables,stddev.participants,r+4)
Tables <- insertrow(Tables,newline,r+5)
Tables <- insertrow(Tables,number.of.scenarios,r+6)
Tables <- insertrow(Tables,stddev.scenarios,r+7)
Tables <- insertrow(Tables, hline, r+8)
# Write the table to a file that can be inserted into the document
write.table(Tables,file=here(".." ,"..", "Manuscript", "Tables", "RETable(Sources).tex"),
            sep="",row.names= FALSE,na="", quote = FALSE, col.names = FALSE)




